module EQSignal

using DSP

export EQSignal

export accInt, acc2vd, baseLineCorr
export targetDispCorr, targetDispVelCorr, targetDispVelAccCorr
export setupSpectra, setSpectraPeriods
export setTargetSpectra, fitSpectra, adjustSpectra, creatArtifitialEQ
export saveEQSignal, saveSpectra

require("EQSignal_h.jl")
require("utility.jl")
# using DataFrames

function Spectra{T<:Real}(SPT::AbstractArray{T,1}, SPA::AbstractArray{T,1}, SPV::AbstractArray{T,1}, SPD::AbstractArray{T,1}, INDSPA::AbstractArray{Integer,1}, SPAT::AbstractArray{T,1},Zeta::T=0.05,SMethod::String="mixed",pseudo::Bool=false)

    NT::Integer = length(SPT)
    Tstart:: T  = SPT[1]
    Tstop :: T  = SPT[end]

    DMode::String = "user"

    return Spectra(Tstart, Tstop, NT, DMode, SPT, SPA, SPV, SPD, INDSPA, SPAT, Zeta, SMethod, pseudo)
end

function Spectra{T<:Real,S<:Integer}(Tstart::T=0.0,Tstop::T=10.0,NT::S=75,DMode::String="log",Zeta::T=0.05,SMethod::String="mixed",pseudo::Bool=false)

    if DMode == "linear"
        SPT = linspace(Tstart,Tstop,NT)
    elseif DMode == "log"
        Tstart = Tstart==0.0 ? 0.02 : Tstart
        SPT = logspace(log10(Tstart),log10(Tstop),NT)
    end

    SPA    = zeros(T,NT)
    SPV    = zeros(T,NT)
    SPD    = zeros(T,NT)
    INDSPA = ones(Integer,NT)
    SPAT   = zeros(T,NT)

    return Spectra(Tstart, Tstop, NT, DMode, SPT, SPA, SPV, SPD, INDSPA, SPAT, Zeta, SMethod, pseudo)
end

function setPeriods{T<:Real}(sp::Spectra,SPT::AbstractArray{T,1})
    NT = length(SPT)
    sp.NT = NT
    sp.DMode = "user"
    sp.SPT = SPT
    sp.Tstart = SPT[1]
    sp.Tstop = SPT[end]
    sp.SPA = zeros(T,NT)
    sp.SPV = zeros(T,NT)
    sp.SPD = zeros(T,NT)
    sp.SPAT = zeros(T,NT)
    nothing
end

function setPeriods(sp::Spectra,NT::Integer)
    if sp.DMode == "linear"
        SPT = linspace(sp.Tstart,sp.Tstop,NT)
    elseif sp.DMode == "log"
        sp.Tstart = sp.Tstart==0.0 ? 0.01 : sp.Tstart
        SPT = logspace(log10(sp.Tstart),log10(sp.Tstop),NT)
    end
    sp.NT = NT
    sp.SPT = SPT
    sp.SPA = zeros(Float64,NT)
    sp.SPV = zeros(Float64,NT)
    sp.SPD = zeros(Float64,NT)
    sp.SPAT = zeros(Float64,NT)
    nothing
end

function setPeriods{T<:Real}(sp::Spectra,NT::Integer,Tstart::T,Tstop::T,DMode::String="log")
    sp.Tstart = Tstart
    sp.Tstop = Tstop
    sp.DMode = DMode

    if sp.DMode == "linear"
        SPT = linspace(sp.Tstart,sp.Tstop,NT)
    elseif sp.DMode == "log"
        sp.Tstart = sp.Tstart==0.0 ? 0.01 : sp.Tstart
        SPT = logspace(log10(sp.Tstart),log10(sp.Tstop),NT)
    end

    sp.NT = NT
    sp.SPT = SPT
    sp.SPA = zeros(T,NT)
    sp.SPV = zeros(T,NT)
    sp.SPD = zeros(T,NT)
    sp.SPAT = zeros(T,NT)

    nothing
end

function EQSignal(name::String,acc::Vector{Float64},t::Vector{Float64},v0::Float64,d0::Float64)

    n = length(acc)

    vel  = zeros(n)
    disp = zeros(n)

    dt = t[2]-t[1]

    Ia    = zeros(n)
    acc0  = deepcopy(acc)
    vel0  = zeros(n)
    disp0 = zeros(n)
    acc_raw  = deepcopy(acc)
    vel_raw  = zeros(n)
    disp_raw = zeros(n)
    paras = Dict{Symbol,Real}()
    for key = [:ARMS,:VRMS,:DRMS]
        paras[key] = 0.0
    end
    sp = Spectra[Spectra()]

    EQSignal(name,acc,vel,disp,t,dt,n,v0,d0,Ia,acc0,vel0,disp0,acc_raw,vel_raw,disp_raw,paras,sp)

end

function EQSignal(name::String,acc::Vector{Float64},dt::Float64,v0::Float64,d0::Float64)

    n = length(acc)

    vel = zeros(n)
    disp = zeros(n)

    t = linspace(0.0,(n-1)*dt,n)

    Ia = zeros(n)
    acc0 = deepcopy(acc)
    vel0 = zeros(n)
    disp0= zeros(n)
    acc_raw  = deepcopy(acc)
    vel_raw  = zeros(n)
    disp_raw = zeros(n)
    paras = Dict{Symbol,Real}()
    for key = [:ARMS,:VRMS,:DRMS]
        paras[key] = 0.0
    end

    sp = Spectra[Spectra()]

    EQSignal(name,acc,vel,disp,t,dt,n,v0,d0,Ia,acc0,vel0,disp0,acc_raw,vel_raw,disp_raw,paras,sp)

end

EQSignal(name::String,acc::Vector{Float64},dt::Float64) = EQSignal(name,acc,dt,0.0,0.0)
EQSignal(name::String,acc::Vector{Float64},t::Vector{Float64}) = EQSignal(name,acc,t,0.0,0.0)
EQSignal(name::String,N::Int64,dt::Float64) = EQSignal(name,zeros(N),dt,0.0,0.0)

function EQSignal(name::String,dt::Float64=0.02)
    
    fn = "$name.txt"
    data = loadtxt(fn)

    if length(size(data)) == 1
        acc = data[:]
        return EQSignal(name,acc,dt)
    else
        t   = data[:,1]
        acc = data[:,2]
        return EQSignal(name,acc,t)
    end

end

EQSignal(acc::Vector{Float64},dt::Float64) = EQSignal("eqs",acc,dt,0.0,0.0)
EQSignal(acc::Vector{Float64},t::Vector{Float64}) = EQSignal("eqs",acc,t,0.0,0.0)
EQSignal(N::Int64,dt::Float64) = EQSignal("eqs",zeros(N),dt,0.0,0.0)
EQSignal(t::Vector{Float64}) = EQSignal("eqs",zeros(length(t)),t[2]-t[1],0.0,0.0)
EQSignal() = EQSignal("eqs",zeros(2048),0.02,0.0,0.0)

function EQSignal(rsn::Integer,dir::Integer=1)
    url_pre = "http://peer.berkeley.edu/nga_files/ath/"
    ngadb,header = readdlm("nga.csv",',';header=true)

    ffn = ngadb[rsn,dir+9]
    fn = download(url_pre*ffn,split(ffn,"/")[2])
    name = split(fn,".")[1]

    acc,N,dt = readNGA(fn)

    return EQSignal(name,acc,dt)
end

# function EQSignal(rsn::Integer,dir::Integer=1)
#     url_pre = "http://peer.berkeley.edu/nga_files/ath/"
#     EQDIR = Symbol[:FNH1,:FNH2,:FNV]
#     ngadb = readtable("nga.csv")

#     ffn = ngadb[EQDIR[dir]][rsn]
#     fn = download(url_pre*ffn,split(ffn,"/")[2])
#     name = split(fn,".")[1]

#     acc,N,dt = readNGA(fn)

#     return EQSignal(name,acc,dt)
# end

function setSpectraPeriods{T<:Real}(eqs::EQSignal,SPT::AbstractArray{T,1})

    for i = 1:length(eqs.sp)
        setPeriods(eqs.sp[i],SPT)
    end
    nothing
end

function setSpectraPeriods(eqs::EQSignal,NT::Integer)

    for i = 1:length(eqs.sp)
        setPeriods(eqs.sp[i],NT)
    end
    nothing
end

function setSpectraPeriods{T<:Real}(eqs::EQSignal,NT::Integer,Tstart::T,Tstop::T,DMode::String="log")

    for i = 1:length(eqs.sp)
        setPeriods(eqs.sp[i],NT,Tstart,Tstop,DMode)
    end
    nothing
end

function setSpectraPara(eqs::EQSignal;Zeta::Float64=0.05,SMethod::String="mixed",pseudo::Bool=false)

    for i = 1:length(eqs.sp)
        eqs.sp[i].zeta = Zeta
        eqs.sp[i].SMethod = SMethod
        eqs.sp[i].pseudo = pseudo
    end
    nothing
end

function setupSpectra{T<:Real}(eqs::EQSignal,Zeta::AbstractArray{T,1};Tstart::Float64=0.02,Tstop::Float64=10.0,NT::Integer=75,DMode::String="log",SMethod::String="mixed",pseudo::Bool=false)

    eqs.sp = Spectra[Spectra(Tstart,Tstop,NT,DMode,zi,SMethod,pseudo) for zi in Zeta]
    nothing
end

function getSpectra(eqs::EQSignal)

    for i = 1:length(eqs.sp)
        eqs.sp[i].SPA,eqs.sp[i].SPV,eqs.sp[i].SPD,eqs.sp[i].INDSPA = responseSpectrum(eqs.acc,eqs.sp[i].SPT,eqs.dt,eqs.sp[i].zeta;method=eqs.sp[i].SMethod,ifpseudo=eqs.sp[i].pseudo,ifabs=true)
    end

end

function getSpectra{T<:Real,S<:Integer}(eqs::EQSignal, SPA::AbstractArray{T,1}, SPV::AbstractArray{T,1}, SPD::AbstractArray{T,1}, INDSPA::AbstractArray{S,1})

    for i = 1:length(eqs.sp)
        eqs.sp[i].SPA,eqs.sp[i].SPV,eqs.sp[i].SPD,eqs.sp[i].INDSPA = abs(SPA), abs(SPV), abs(SPD), INDSPA
    end
    nothing
end

function setTargetSpectra{T<:Real}(eqs::EQSignal,Tg::T=0.9,amax::T=2.25;code="GB50011-2010")

    for i = 1:length(eqs.sp)
        SPAT, = codeSpectrum(Tg,eqs.sp[i].SPT,eqs.sp[i].zeta;amax=amax)
        eqs.sp[i].SPAT = SPAT
    end
    nothing
end

function getEPA{T<:Real}(eqs::EQSignal,P0::T=0.1,P1::T=0.9,AMAX::T=2.25,spid::Integer=1)
    INDP = find(p->p>=P0&&p<=P1,eqs.sp.SPT)
    return trapz(eqs.sp[spid].SPT[INDP],eqs.sp[spid].SPA[INDP])/(P1-P0)/AMAX
end

function fitSpectra(eqs::EQSignal,spid::Integer=1;tol::Float64=0.05,maxiter::Integer=15)

    sp = eqs.sp[spid]
    acc,SPA,SPV,SPD,INDSPA = SPFit(eqs.acc,sp.SPT,sp.SPAT,eqs.dt,sp.zeta;method=sp.SMethod,tol=tol,maxiter=maxiter)
    eqs.acc = acc
    eqs.acc0 = deepcopy(acc)
    getSpectra(eqs,SPA,SPV,SPD,INDSPA)
    accInt(eqs)
    nothing
end

function adjustSpectra(eqs::EQSignal,spid::Integer=1;tol::Float64=0.05,maxiter::Integer=5)
    sp = eqs.sp[spid]
    # getSpectra(eqs)
    acc,SPA,SPV,SPD,INDSPA = SPAdjust(eqs.acc,sp.SPT,sp.SPAT,sp.zeta,eqs.dt;method=sp.SMethod,tol=tol,maxiter=maxiter)
    eqs.acc = acc
    eqs.acc0 = deepcopy(acc)
    getSpectra(eqs,abs(SPA),abs(SPV),abs(SPD),INDSPA)
    accInt(eqs)
    nothing
end

function adjustSpectra(eqs::EQSignal,spid::AbstractArray{Integer,1};tol::Float64=0.05,maxiter::Integer=5)
    sps = eqs.sp[spid]
    Zeta = Float64[spi.zeta for spi in sps]
    # getSpectra(eqs)
    SPAT, = codeSpectrum(Tg,sps[1].SPT,Zeta)
    acc,SPA,SPV,SPD,INDSPA = SPAdjust(eqs.acc,sps[1].SPT,SPAT,Zeta,eqs.dt;method=sp.SMethod,tol=tol,maxiter=maxiter)
    eqs.acc = acc
    eqs.acc0 = deepcopy(acc)
    # getSpectra(eqs,abs(SPA),abs(SPV),abs(SPD),INDSPA)
    accInt(eqs)
    nothing
end

function EQSignal()
    
end

function creatArtifitialEQ(eqs::EQSignal,spid::Integer=1;tol::Float64=0.05,maxiter::Integer=5)
    evlp = Envelope(eqs.t)
    sp = eqs.sp[spid]
    acc,SPA,SPV,SPD,INDSPA = ArtificialWave(eqs.t,evlp,sp.SPT,sp.SPAT,sp.zeta;ifdiffphase=false,method=sp.SMethod,tol=tol,maxiter=maxiter)
    eqs.acc = acc
    eqs.acc0 = deepcopy(acc)
    getSpectra(eqs,SPA,SPV,SPD,INDSPA)
    accInt(eqs)
    nothing
end

# function initialize(eqs::EQSignal)
#     setInitial(eqs)
#     accInt(eqs)
#     rationalAccInt(eqs)
#     getRMS(eqs)
#     setInitial(eqs)
# end

function Normalize(eqs::EQSignal)
    eqs.acc = eqs.acc ./ peak(eqs.acc)
    nothing
end

function accInt(eqs::EQSignal;original::Bool=false)

    if original
        eqs.vel  = cumtrapz(eqs.t,eqs.acc0,eqs.v0)
    else
        eqs.vel  = cumtrapz(eqs.t,eqs.acc,eqs.v0)
    end

    eqs.disp = cumtrapz(eqs.t,eqs.vel,eqs.d0)
    nothing

end

acc2vd(eqs::EQSignal) = accInt(eqs)

function adjustAccPeak(eqs::EQSignal,pv::Float64=1.0)
    adjustpeak!(eqs.acc)
end

function setCurrentAsOriginal(eqs::EQSignal)
    eqs.acc0  = deepcopy(eqs.acc)
    eqs.vel0  = deepcopy(eqs.vel)
    eqs.disp0 = deepcopy(eqs.disp)
end

function setCurrentAsRaw(eqs::EQSignal)
    eqs.acc_raw  = deepcopy(eqs.acc)
    eqs.vel_raw  = deepcopy(eqs.vel)
    eqs.disp_raw = deepcopy(eqs.disp)
end

function getRMS(eqs::EQSignal)

    eqs.paras[:ARMS] = rms(eqs.acc0)
    eqs.paras[:VRMS] = rms(eqs.vel0)
    eqs.paras[:DRMS] = rms(eqs.disp0)

    nothing

end

function getAriasIntensity(eqs::EQSignal)

    eqs.Ia = cumtrapz(eqs.t,0.5*pi*eqs.acc.^2)
    nothing

end

function getSignificantDuration(eqs::EQSignal,thresholds::Vector{Float64}=[0.05,0.95];output::String="time")

    IA = eqs.Ia[end]
    tstart,tend = 0.0,eqs.t[end]
    idxstart,idxend = 1,eqs.n

    for i = 1:eqs.n
        if eqs.Ia[i] > thresholds[1]*IA
            tstart   = eqs.t[i]
            idxstart = i
            break
        end
    end

    for i = 1:eqs.n
        if eqs.Ia[end-i+1] < thresholds[2]*IA
            tend   = eqs.t[end-i+1]
            idxend = eqs.n-i+1
            break
        end
    end

    if output == "index"
        return idxstart,idxend
    elseif output == "time"
        return tstart,tend
    end
end

function trim(eqs::EQSignal,idxstart::Int64,idxend::Int64)

    eqs.n = idxend-idxstart+1
    eqs.t = linspace(0.0,(eqs.n-1)*eqs.dt,eqs.n)

    eqs.acc   = eqs.acc[idxstart:idxend]
    eqs.vel   = eqs.vel[idxstart:idxend]
    eqs.disp  = eqs.disp[idxstart:idxend]
    eqs.acc0  = eqs.acc0[idxstart:idxend]
    eqs.vel0  = eqs.vel0[idxstart:idxend]
    eqs.disp0 = eqs.disp0[idxstart:idxend]

    eqs.Ia    = zeros(eqs.n)
    # accInt(eqs)

end

function recover(eqs::EQSignal)

    eqs.acc = deepcopy(eqs.acc0)
    setInitial(eqs)

end

function setInitial(eqs::EQSignal,v0::Float64=0.0,d0::Float64=0.0)

    eqs.v0 = v0
    eqs.d0 = d0

    nothing
end

function getRationalInitial(eqs::EQSignal,tk::Float64,vk::Float64,dk::Float64)

    ik = int(tk/eqs.dt)
    vel  = cumtrapz(eqs.t[1:ik],eqs.acc0[1:ik],eqs.v0)
    disp  = cumtrapz(eqs.t[1:ik],vel,eqs.d0)

    eqs.v0 = vk - vel[end]
    eqs.d0 = dk - eqs.v0*(tk-eqs.t[1]) - disp[end]

end

function getRationalInitial(eqs::EQSignal)

    accInt(eqs;original=true)

    cv0 = mean(eqs.vel)
    cd0 = mean(eqs.disp .- cv0*eqs.t)
    eqs.v0 = eqs.v0 - cv0
    eqs.d0 = eqs.d0 - cd0

end

function getTargetDispVelAcc(eqs::EQSignal)

    accInt(eqs;original=true)

    # aslope = (eqs.acc[end]-0.0)/(eqs.t[end]-eqs.t[1])
    # eqs.acc0 = eqs.acc .- aslope.*(eqs.t.-eqs.t[1])

    vslope = (eqs.vel[end]-0.0)/(eqs.t[end]-eqs.t[1])
    eqs.vel0 = eqs.vel .- vslope.*(eqs.t.-eqs.t[1])

    dslope = (eqs.disp[end]-0.0)/(eqs.t[end]-eqs.t[1])
    eqs.disp0 = eqs.disp .- dslope.*(eqs.t.-eqs.t[1])
    nothing

end

function rationalAccInt(eqs::EQSignal)
    # calculate the "real" initial vel & disp and then integrate acc using the "real" initial value
    v0, d0 = eqs.v0, eqs,d0
    getRationalInitial(eqs)
    eqs.vel0 = cumtrapz(eqs.t,eqs.acc0,eqs.v0)
    eqs.disp0= cumtrapz(eqs.t,eqs.vel0,eqs.d0)
    setInitial(eqs,v0,d0)
    nothing
end

function detrendDisp(eqs::EQSignal,oh::Int64=1,ol::Int64=1)
    eqs.disp = detrend(eqs.disp,oh,ol)
    nothing
end

function detrendVel(eqs::EQSignal,oh::Int64=0,ol::Int64=0)
    eqs.vel = detrend(eqs.vel,oh,ol)
    nothing
end

function detrendAcc(eqs::EQSignal,oh::Int64=3,ol::Int64=0)
    eqs.acc = detrend(eqs.acc,oh,ol)
    accInt(eqs)
    nothing
end

function stepDetrendAcc(eqs::EQSignal,ns::Int64=4,oh::Int64=3,ol::Int64=0)
    acc = eqs.acc
    sidx = int(linspace(eqs.n/ns,eqs.n,ns))

    for s = sidx
        acc[1:s] = detrend(acc[1:s],oh,ol)
    end
    accInt(eqs)
    nothing
end

function baseLineCorr(eqs::EQSignal,oh::Int64=4,ol::Int64=2)

    accInt(eqs;original=true)
    eqs.acc = basecorrv(eqs.acc,eqs.vel,eqs.t,max(oh-1,0),max(ol-1,0))
    accInt(eqs)
    eqs.acc = basecorrd(eqs.acc,eqs.disp,eqs.t,oh,ol)
    accInt(eqs)

end

function stepBaseLineCorr(eqs::EQSignal,ns::Int64=4,oh::Int64=4,ol::Int64=2)

    accInt(eqs;original=true)
    if ns == 1
        baseLineCorr(eqs,oh,ol)
        return nothing
    else
        sidx = round(Int64,linspace(eqs.n/ns,eqs.n,ns))
    end

    for s = sidx
        eqs.acc[1:s] = basecorrv(eqs.acc[1:s],eqs.vel[1:s],eqs.t[1:s],max(oh-1,0),max(ol-1,0))
        accInt(eqs)
        eqs.acc[1:s] = basecorrd(eqs.acc[1:s],eqs.disp[1:s],eqs.t[1:s],oh,ol)
        accInt(eqs)
    end
end

function DriftRatio(eqs::EQSignal)
    rd = eqs.paras[:DRMS]
    md = mean(eqs.disp)
    abs(md/rd)
end

function AmpRatio(eqs::EQSignal)
    rd = eqs.paras[:DRMS]
    ad = rms(eqs.disp)
    abs(ad/rd)
end

function highpassFilt(eqs::EQSignal)

end

function finalVelZero(eqs::EQSignal)
    eqs.acc = eqs.acc .- eqs.vel[end]/eqs.t[end]
    accInt(eqs)
end

function finalDispZero(eqs::EQSignal)
    eqs.acc = eqs.acc .- 2.0*eqs.disp[end]/eqs.t[end]^2
    accInt(eqs)
end

function finalVelDispZero(eqs::EQSignal)
    const0  = -6.0*eqs.disp[end]/eqs.t[end]^2+2.0*eqs.vel[end]/eqs.t[end]
    linear  = 12.0*eqs.disp[end]/eqs.t[end]^3-6.0*eqs.vel[end]/eqs.t[end]^2
    corr    = const0 .+ linear*eqs.t
    eqs.acc = eqs.acc .+ corr
    accInt(eqs)
end

function targetDispCorr(eqs::EQSignal,dispt::Vector{Float64},np::Int64=4,nm::Int64=4)

    if nm > np
        nm = np
    end

    if np == 1
        cidx = int([eqs.n*0.618])
        # cidx = round(Int,[eqs.n*0.618])
    else
        cidx = int(linspace(eqs.n/np,eqs.n,np))
        # cidx = round(Int,linspace(eqs.n/np,eqs.n,np))
    end

    # f(i::Int64,t::Float64)   = t^(i-1)
    # fp(i::Int64,t::Float64)  = i>1 ? (i-1)*t^(i-2): 0.0
    # fpp(i::Int64,t::Float64) = i>2 ? (i-1)*(i-2)*t^(i-3): 0.0

    f(i::Int64,t::Float64)   = t^(i+1)
    fp(i::Int64,t::Float64)  = i>-1 ? (i+1)*t^(i): 0.0
    fpp(i::Int64,t::Float64) = i>0  ? (i+1)*(i)*t^(i-1): 0.0

    accInt(eqs;original=true)
    disp0 = deepcopy(eqs.disp) # save original disp

    tc = zeros(np) # controlled time
    b = zeros(np) # right-hand-side vector
    A = zeros(np,nm) # coeff matrix

    # setup tc, b, A
    for k = 1:np
        tc[k] = eqs.t[cidx[k]]
        b[k]  = dispt[cidx[k]]-disp0[cidx[k]]
        for i = 1:nm
            A[k,i] = f(i,tc[k])
        end
    end

    # solve linear system to get c[i]
    c = A\b

    corr = zeros(eqs.n) # correction time series

    # calculate corr using solved c[i]
    for k = 1:eqs.n
        tk = eqs.t[k]
        for i = 1:nm
            corr[k] += c[i]*fpp(i,tk)
        end
    end

    eqs.acc = eqs.acc .+ corr
    accInt(eqs)
    nothing

end

function targetDispCorr(eqs::EQSignal,dispt::Vector{Float64},vend::Float64,np::Int64=4,nm::Int64=4)

    if nm > np
        nm = np
    end

    nm += 1

    if np == 1
        cidx = int([eqs.n*0.618])
        # cidx = round(Int,[eqs.n*0.618])
    else
        cidx = int(linspace(eqs.n/np,eqs.n,np))
        # cidx = round(Int,linspace(eqs.n/np,eqs.n,np))
    end

    # f(i::Int64,t::Float64)   = t^(i-1)
    # fp(i::Int64,t::Float64)  = i>1 ? (i-1)*t^(i-2): 0.0
    # fpp(i::Int64,t::Float64) = i>2 ? (i-1)*(i-2)*t^(i-3): 0.0

    f(i::Int64,t::Float64)   = t^(i+1)
    fp(i::Int64,t::Float64)  = i>-1 ? (i+1)*t^(i): 0.0
    fpp(i::Int64,t::Float64) = i>0  ? (i+1)*(i)*t^(i-1): 0.0

    accInt(eqs;original=true)
    disp0 = deepcopy(eqs.disp) # save original disp
    vel0  = deepcopy(eqs.vel)  # save original vel

    tc = zeros(np+1) # controlled time
    b = zeros(np+1) # right-hand-side vector
    A = zeros(np+1,nm) # coeff matrix

    # setup tc, b, A
    for k = 1:np
        tc[k] = eqs.t[cidx[k]]
        b[k]  = dispt[cidx[k]]-disp0[cidx[k]]
        for i = 1:nm
            A[k,i] = f(i,tc[k])
        end
    end

    # final vel constraint
    tc[np+1] = eqs.t[end]
    b[np+1] = vend-vel0[end]
    for i = 1:nm
        A[np+1,i] = fp(i,tc[np+1])
    end

    # solve linear system to get c[i]
    c = A\b

    # initial correction time series
    corr = zeros(eqs.n)

    # calculate corr using solved c[i]
    for k = 1:eqs.n
        tk = eqs.t[k]
        for i = 1:nm
            corr[k] += c[i]*fpp(i,tk)
        end
    end

    # apply correction
    eqs.acc = eqs.acc .+ corr
    accInt(eqs)
    nothing

end

function targetDispCorr(eqs::EQSignal,dispt::Vector{Float64},vend::Float64,aend::Float64,np::Int64=4,nm::Int64=4)

    if nm > np
        nm = np
    end

    nm += 2

    if np == 1
        cidx = int([eqs.n*0.618])
        # cidx = round(Int,[eqs.n*0.618])
    else
        cidx = int(linspace(eqs.n/np,eqs.n,np))
        # cidx = round(Int,linspace(eqs.n/np,eqs.n,np))
    end

    # f(i::Int64,t::Float64)   = t^(i-1)
    # fp(i::Int64,t::Float64)  = i>1 ? (i-1)*t^(i-2): 0.0
    # fpp(i::Int64,t::Float64) = i>2 ? (i-1)*(i-2)*t^(i-3): 0.0

    f(i::Int64,t::Float64)   = t^(i+1)
    fp(i::Int64,t::Float64)  = i>-1 ? (i+1)*t^(i): 0.0
    fpp(i::Int64,t::Float64) = i>0  ? (i+1)*(i)*t^(i-1): 0.0

    accInt(eqs;original=true)
    disp0 = deepcopy(eqs.disp) # save original disp
    vel0  = deepcopy(eqs.vel)  # save original vel
    acc0  = deepcopy(eqs.acc)  # save original vel

    tc = zeros(np+2) # controlled time
    b = zeros(np+2) # right-hand-side vector
    A = zeros(np+2,nm) # coeff matrix

    # setup tc, b, A
    for k = 1:np
        tc[k] = eqs.t[cidx[k]]
        b[k]  = dispt[cidx[k]]-disp0[cidx[k]]
        for i = 1:nm
            A[k,i] = f(i,tc[k])
        end
    end

    # final vel constraint
    tc[np+1] = eqs.t[end]
    b[np+1] = vend-vel0[end]
    for i = 1:nm
        A[np+1,i] = fp(i,tc[np+1])
    end
    # final acc constraint
    tc[np+2] = eqs.t[end]
    b[np+2] = aend-acc0[end]
    for i = 1:nm
        A[np+2,i] = fpp(i,tc[np+2])
    end

    # solve linear system to get c[i]
    c = A\b

    # initial correction time series
    corr = zeros(eqs.n)

    # calculate corr using solved c[i]
    for k = 1:eqs.n
        tk = eqs.t[k]
        for i = 1:nm
            corr[k] += c[i]*fpp(i,tk)
        end
    end

    # apply correction
    eqs.acc = eqs.acc .+ corr
    accInt(eqs)
    nothing

end

function targetDispCorr(eqs::EQSignal,np::Int64=4,nm::Int64=4)

    # getRationalInitial(eqs)
    # accInt(eqs)
    # setInitial(eqs)

    dispt = eqs.disp0
    targetDispCorr(eqs,dispt,np,nm)

end

function targetDispCorr(eqs::EQSignal,vend::Float64,np::Int64=4,nm::Int64=4)

    dispt = eqs.disp0
    targetDispCorr(eqs,dispt,vend,np,nm)

end

function targetDispCorr(eqs::EQSignal,vend::Float64,aend::Float64,np::Int64=4,nm::Int64=4)

    dispt = eqs.disp0
    targetDispCorr(eqs,dispt,vend,aend,np,nm)
end

function targetDispVelCorr(eqs::EQSignal,dispt::Vector{Float64},velt::Vector{Float64},np::Int64=4,nm::Int64=8)

    if nm > 2*np
        nm = 2*np
    end

    if np == 1
        cidx = int([eqs.n*1.0])
        # cidx = round(Int,[eqs.n*0.618])
    else
        cidx = int(linspace(eqs.n/np,eqs.n,np))
        # cidx = round(Int,linspace(eqs.n/np,eqs.n,np))
    end

    # f(i::Int64,t::Float64)   = t^(i-1)
    # fp(i::Int64,t::Float64)  = i>1 ? (i-1)*t^(i-2): 0.0
    # fpp(i::Int64,t::Float64) = i>2 ? (i-1)*(i-2)*t^(i-3): 0.0

    # f(i::Int64,t::Float64)   = t^(i+1)
    # fp(i::Int64,t::Float64)  = i>-1 ? (i+1)*t^(i): 0.0
    # fpp(i::Int64,t::Float64) = i>0  ? (i+1)*(i)*t^(i-1): 0.0

    f(i::Int64,t::Float64)   = t^(i+1+1)
    fp(i::Int64,t::Float64)  = i+1>-1 ? (i+1+1)*t^(i+1): 0.0
    fpp(i::Int64,t::Float64) = i+1>0  ? (i+1+1)*(i+1)*t^(i+1-1): 0.0

    accInt(eqs;original=true)
    disp0 = deepcopy(eqs.disp) # save original disp
    vel0  = deepcopy(eqs.vel)  # save original vel
    acc0  = deepcopy(eqs.acc)  # save original vel

    tc = zeros(np) # controlled time
    b = zeros(np*2) # right-hand-side vector
    A = zeros(np*2,nm) # coeff matrix

    # setup tc, b, A
    for k = 1:np
        tc[k] = eqs.t[cidx[k]]
        b[k]  = dispt[cidx[k]]-disp0[cidx[k]]
        b[k+np]  = velt[cidx[k]]-vel0[cidx[k]]
        for i = 1:nm
            A[k,i] = f(i,tc[k])
            A[k+np,i] = fp(i,tc[k])
        end
    end

    # # final vel constraint
    # tc[np+1] = eqs.t[end]
    # b[np+1] = vend-vel0[end]
    # for i = 1:nm
    #     A[np+1,i] = fp(i,tc[np+1])
    # end
    # # final acc constraint
    # tc[np+2] = eqs.t[end]
    # b[np+2] = aend-acc0[end]
    # for i = 1:nm
    #     A[np+2,i] = fpp(i,tc[np+2])
    # end

    # solve linear system to get c[i]
    c = A\b

    # initial correction time series
    corr = zeros(eqs.n)

    # calculate corr using solved c[i]
    for k = 1:eqs.n
        tk = eqs.t[k]
        for i = 1:nm
            corr[k] += c[i]*fpp(i,tk)
        end
    end

    # apply correction
    eqs.acc = eqs.acc .+ corr
    accInt(eqs)
    nothing

end

function targetDispVelCorr(eqs::EQSignal,np::Int64=2,nm::Int64=4)

    dispt = eqs.disp0
    velt = eqs.vel0
    targetDispVelCorr(eqs,dispt,velt,np,nm)

end

function targetDispVelCorr(eqs::EQSignal,dispt::Vector{Float64},velt::Vector{Float64},aend::Float64,np::Int64=2,nm::Int64=4)

    if nm > 2*np
        nm = 2*np
    end

    nm += 1

    if np == 1
        cidx = int([eqs.n*1.0])
        # cidx = round(Int,[eqs.n*0.618])
    else
        cidx = int(linspace(eqs.n/np,eqs.n,np))
        # cidx = round(Int,linspace(eqs.n/np,eqs.n,np))
    end

    # f(i::Int64,t::Float64)   = t^(i-1)
    # fp(i::Int64,t::Float64)  = i>1 ? (i-1)*t^(i-2): 0.0
    # fpp(i::Int64,t::Float64) = i>2 ? (i-1)*(i-2)*t^(i-3): 0.0

    # f(i::Int64,t::Float64)   = t^(i+1)
    # fp(i::Int64,t::Float64)  = i>-1 ? (i+1)*t^(i): 0.0
    # fpp(i::Int64,t::Float64) = i>0  ? (i+1)*(i)*t^(i-1): 0.0

    f(i::Int64,t::Float64)   = t^(i+1+1)
    fp(i::Int64,t::Float64)  = i+1>-1 ? (i+1+1)*t^(i+1): 0.0
    fpp(i::Int64,t::Float64) = i+1>0  ? (i+1+1)*(i+1)*t^(i+1-1): 0.0

    accInt(eqs;original=true)
    disp0 = deepcopy(eqs.disp) # save original disp
    vel0  = deepcopy(eqs.vel)  # save original vel
    acc0  = deepcopy(eqs.acc)  # save original vel

    tc = zeros(np+1) # controlled time
    b = zeros(np*2+1) # right-hand-side vector
    A = zeros(np*2+1,nm) # coeff matrix

    # setup tc, b, A
    for k = 1:np
        tc[k] = eqs.t[cidx[k]]
        b[k]  = dispt[cidx[k]]-disp0[cidx[k]]
        b[k+np]  = velt[cidx[k]]-vel0[cidx[k]]
        for i = 1:nm
            A[k,i] = f(i,tc[k])
            A[k+np,i] = fp(i,tc[k])
        end
    end

    # # final vel constraint
    # tc[np+1] = eqs.t[end]
    # b[np+1] = vend-vel0[end]
    # for i = 1:nm
    #     A[np+1,i] = fp(i,tc[np+1])
    # end
    # final acc constraint
    tc[np+1] = eqs.t[end]
    b[np*2+1] = aend-acc0[end]
    for i = 1:nm
        A[np*2+1,i] = fpp(i,tc[np+1])
    end

    # solve linear system to get c[i]
    c = A\b

    # initial correction time series
    corr = zeros(eqs.n)

    # calculate corr using solved c[i]
    for k = 1:eqs.n
        tk = eqs.t[k]
        for i = 1:nm
            corr[k] += c[i]*fpp(i,tk)
        end
    end

    # apply correction
    eqs.acc = eqs.acc .+ corr
    accInt(eqs)
    nothing

end

function targetDispVelCorr(eqs::EQSignal,aend::Float64,np::Int64=2,nm::Int64=4)

    dispt = eqs.disp0
    velt = eqs.vel0
    targetDispVelCorr(eqs,dispt,velt,aend,np,nm)
end

function targetDispVelAccCorr(eqs::EQSignal,dispt::Vector{Float64},velt::Vector{Float64},acct::Vector{Float64},np::Int64=2,nm::Int64=6)

    if nm > 3*np
        nm = 3*np
    end

    if np == 1
        cidx = int([eqs.n*0.618])
        # cidx = round(Int,[eqs.n*0.618])
    else
        cidx = int(linspace(eqs.n/np,eqs.n,np))
        # cidx = round(Int,linspace(eqs.n/np,eqs.n,np))
    end

    # f(i::Int64,t::Float64)   = t^(i-1)
    # fp(i::Int64,t::Float64)  = i>1 ? (i-1)*t^(i-2): 0.0
    # fpp(i::Int64,t::Float64) = i>2 ? (i-1)*(i-2)*t^(i-3): 0.0

    # f(i::Int64,t::Float64)   = t^(i+1)
    # fp(i::Int64,t::Float64)  = i>-1 ? (i+1)*t^(i): 0.0
    # fpp(i::Int64,t::Float64) = i>0  ? (i+1)*(i)*t^(i-1): 0.0

    f(i::Int64,t::Float64)   = t^(i+1+1)
    fp(i::Int64,t::Float64)  = i+1>-1 ? (i+1+1)*t^(i+1): 0.0
    fpp(i::Int64,t::Float64) = i+1>0  ? (i+1+1)*(i+1)*t^(i+1-1): 0.0

    accInt(eqs;original=true)
    disp0 = deepcopy(eqs.disp) # save original disp
    vel0  = deepcopy(eqs.vel)  # save original vel
    acc0  = deepcopy(eqs.acc)  # save original vel

    tc = zeros(np) # controlled time
    b = zeros(np*3) # right-hand-side vector
    A = zeros(np*3,nm) # coeff matrix

    # setup tc, b, A
    for k = 1:np
        tc[k] = eqs.t[cidx[k]]
        b[k]  = dispt[cidx[k]]-disp0[cidx[k]]
        b[k+np]  = velt[cidx[k]]-vel0[cidx[k]]
        b[k+2*np]  = acct[cidx[k]]-acc0[cidx[k]]
        for i = 1:nm
            A[k,i] = f(i,tc[k])
            A[k+np,i] = fp(i,tc[k])
            A[k+2*np,i] = fpp(i,tc[k])
        end
    end

    # # final vel constraint
    # tc[np+1] = eqs.t[end]
    # b[np+1] = vend-vel0[end]
    # for i = 1:nm
    #     A[np+1,i] = fp(i,tc[np+1])
    # end
    # # final acc constraint
    # tc[np+2] = eqs.t[end]
    # b[np+2] = aend-acc0[end]
    # for i = 1:nm
    #     A[np+2,i] = fpp(i,tc[np+2])
    # end

    # solve linear system to get c[i]
    c = A\b

    # initial correction time series
    corr = zeros(eqs.n)

    # calculate corr using solved c[i]
    for k = 1:eqs.n
        tk = eqs.t[k]
        for i = 1:nm
            corr[k] += c[i]*fpp(i,tk)
        end
    end

    # apply correction
    eqs.acc = eqs.acc .+ corr
    accInt(eqs)
    nothing

end

function targetDispVelAccCorr(eqs::EQSignal,np::Int64=2,nm::Int64=6)

    dispt = eqs.disp0
    velt = eqs.vel0
    acct = eqs.acc0
    targetDispVelAccCorr(eqs,dispt,velt,acct,np,nm)

end

function beginWithZero(eqs::EQSignal, Nc::Integer)

    a0 = eqs.acc[1]
    aslope = (0.0-a0)/(eqs.t[Nc]-eqs.t[1])

    eqs.acc[1:Nc] = eqs.acc[1:Nc] .- aslope.*(eqs.t[1:Nc] .- eqs.t[1]) .- a0

    accInt(eqs)
    nothing
end

function beginWithZero(eqs::EQSignal)

    # Nc = findfirstpeak(eqs.vel)
    Nc = int(0.05*eqs.n)
    beginWithZero(eqs,Nc)

end

function endWithZero(eqs::EQSignal;withAcc::Bool=true)

    # Nc = min(maximum(eqs.sp.INDSPA)+1,int(0.96*eqs.n))
    # Nc = max(Nc,int(0.92*eqs.n))
    Nc = int(0.95*eqs.n)
    endWithZero(eqs,Nc;withAcc=withAcc)
end

function endWithZero(eqs::EQSignal,Nc::Integer;withAcc::Bool=true)

    f(i::Int64,t::Float64)   = t^(i+1+1)
    fp(i::Int64,t::Float64)  = i+1>-1 ? (i+1+1)*t^(i+1): 0.0
    fpp(i::Int64,t::Float64) = i+1>0  ? (i+1+1)*(i+1)*t^(i+1-1): 0.0

    dim = withAcc ? 3 : 2

    tc = eqs.t[end]-eqs.t[Nc] # controlled time
    b = zeros(dim) # right-hand-side vector
    A = zeros(dim,dim) # coeff matrix

    # setup tc, b, A
    b[1] = 0.0 - eqs.disp[end]
    b[2] = 0.0 - eqs.vel[end]
    if withAcc b[3] = 0.0 - eqs.acc[end] end

    for i = 1:dim
        A[1,i] = f(i,tc)
        A[2,i] = fp(i,tc)
        if withAcc A[3,i] = fpp(i,tc) end
    end

    # solve linear system to get c[i]
    c = A\b
    # println((A*c.-b))

    # initial correction time series
    corr = zeros(eqs.n)

    # calculate corr using solved c[i]
    for k = Nc:eqs.n
        tk = eqs.t[k] - eqs.t[Nc]
        for i = 1:dim
            corr[k] += c[i]*fpp(i,tk)
        end
    end

    # apply correction
    eqs.acc = eqs.acc .+ corr
    accInt(eqs)
    nothing
end

function saveEQSignal(eqs::EQSignal)

    output = [eqs.t eqs.acc eqs.vel eqs.disp eqs.Ia/eqs.Ia[end]]
    writedlm("$(eqs.name)-TimeHistory.dat",output)
    # conf = open("fname","a")
    # write(conf,"$(fn)\n")
    # close(conf)
end

function saveSpectra(eqs::EQSignal,fn::String)
    for sp = eqs.sp
        output = [sp.SPT sp.SPA sp.SPV sp.SPD]
        writedlm("$(eqs.name)-%$(int(sp.zeta*100))DampingSpectra.dat",output)
    end
end

# function spectrum(eqs::EQSignal,Periods::AbstractArray{Float64,1};xi::Float64=0.05)
#     SPA,SPV,SPD = similar(Periods),similar(Periods),similar(Periods)
#     for (i,p) = enumerate(Periods)
#         RA,RV,RD,RP,RT = newmarkBeta(eqs.acc,eqs.dt,p,xi)
#         SPA[i],SPV[i],SPD[i] = abs(RP)
#     end
#     return SPA,SPV,SPD
# end

# function spectrum(eqs::EQSignal,NP::Int64;xi::Float64=0.05)
#     Periods = logspace(-2,1,NP)
#     SPA,SPV,SPD = spectrum(eqs,Periods;xi=xi)
#     return Periods,SPA,SPV,SPD
# end

end